Leaflet: An Introduction to Choropleth maps in Rleaflet

Introduction

Leaflet is the leading open-source JavaScript library for mobile-friendly interactive maps. Weighing just about 38 KB of JS, it has all the mapping features most developers ever need.

Leaflet is designed with simplicity, performance and usability in mind. It works efficiently across all major desktop and mobile platforms, can be extended with lots of plugins, has a beautiful, easy to use and well-documented API and a simple, readable source code that is a joy to contribute to.

Leaflet is one of the most popular open-source JavaScript libraries for interactive maps. Its used by websites ranging from The New York Times and The Washington Post to GitHub and Flickr, as well as GIS specialists like OpenStreetMap, Mapbox, and CartoDB.

This R package makes it easy to integrate and control Leaflet maps in R.

A cheatsheet for this package is here: https://ugoproto.github.io/ugo_r_doc/leaflet-cheat-sheet.pdf

#Example of a simple map

#Check if leaflet is installed, if not install it
require("leaflet")

## Loading required package: leaflet

require("DT")

## Loading required package: DT

# Add leaflet extras
require("leaflet.extras")

## Loading required package: leaflet.extras

# A simple map...

m <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=174.768, lat=-36.852, popup="The birthplace of R")
m  # Print the map

Guidance

You can get more information on the using the Leaflet library in r from the Leaflet Github Repository.
You can get more information on Leaflet from the Leaflet website. This github also has some helpful examples

---

Choropleth Map

A choropleth map allows us to view data by shading areas of a map according to values within those areas.

In order to do this we are going to need to add polygons to the map, ready for us to overlay our points. Shading can be applied depending on how many points fall within each polygon.

To demonstrate this, a file has been added to the repo contianing regions of the UK - this is open source and freely available from the ONS.

Getting Started

The first thing to do when creating a choropleth leaflet map using RLeaflet is to load the required packages;

require("readr")

## Loading required package: readr

require("rgdal")

## Loading required package: rgdal

## Loading required package: sp

## rgdal: version: 1.4-3, (SVN revision 828)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Anthony/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Anthony/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.3-1

require("sf")

## Loading required package: sf

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

require("dplyr")

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Data to be plotted

Next, lets read in the accidents data using readr function ‘read_csv’; uk_accidents_2017 <- read_csv(file = “uk_accidents_2017.csv”)

uk_accidents_2017 <- read_csv(file = "uk_accidents_2017.csv")

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   accdate = col_date(format = ""),
##   acctime = col_time(format = ""),
##   rgn = col_character(),
##   lsoa11 = col_character(),
##   msoa11 = col_character(),
##   lsoa_of_acc_loc = col_character()
## )

## See spec(...) for full column specifications.

It’s always a good idea to become familiar with the data;

head(uk_accidents_2017)

## # A tibble: 6 x 26
##   location_eastin~ location_northi~ longitude latitude police_force
##              <dbl>            <dbl>     <dbl>    <dbl>        <dbl>
## 1            91929            12107     -6.29     49.9           50
## 2           135295            26294     -5.70     50.1           50
## 3           136729            26553     -5.68     50.1           50
## 4           136963            31354     -5.68     50.1           50
## 5           137000            30348     -5.68     50.1           50
## 6           137003            31771     -5.68     50.1           50
## # ... with 21 more variables: accident_severity <dbl>,
## #   number_of_vehicles <dbl>, number_of_casualties <dbl>, accdate <date>,
## #   day_of_week <dbl>, acctime <drtn>, rgn <chr>, lsoa11 <chr>,
## #   msoa11 <chr>, road_type <dbl>, speed_limit <dbl>,
## #   junction_detail <dbl>, junction_control <dbl>, light_conditions <dbl>,
## #   weather_conditions <dbl>, road_surface_conditions <dbl>,
## #   special_conditions_at_site <dbl>, carriageway_hazards <dbl>,
## #   urban_or_rural_area <dbl>, did_police_attend <dbl>,
## #   lsoa_of_acc_loc <chr>

Each row in the dataset refers to a single accident that occured somewhere in England. This dataset has some interesting items in it that you may wish to perform analysis on such as; information on severity, date, time, location, number casualties and number of vehicles. There is also an indentifier that can help us link to several of the Ordanance Survey boundaries that you may wish to bring in later.

Lets subset the accidents data for variables we want to visualise on a choropleth map. For this example we’ll just keep our joining variable, ‘rgn’ which links to a region shapefile, and we’ll create a simple count of the number of accidents in each region;

uk_accidents_2017 <- uk_accidents_2017 %>%
  select(rgn) %>%
  group_by(rgn) %>%
  mutate(count = n()) %>%
  distinct()
head(uk_accidents_2017)

## # A tibble: 6 x 2
## # Groups:   rgn [6]
##   rgn       count
##   <chr>     <int>
## 1 E12000009 10055
## 2 W99999999  4554
## 3 E12000002 13164
## 4 E12000005 10977
## 5 E12000001  4088
## 6 E12000003 11336

NOTE: The joining keys ‘rgn’, ‘lsoa11’ and ‘msoa11’ were added to the data prior to this example being compiled.

To recap, we’ve brought in a dataset that contains data we may wish to plot on a map. The data could be used to create points as they have lat/long coordinates but it does not contain the necessary data to plot polygons.

Polygons

Polygons are more complex and require several elements for RLeaflet to be able to plot them, such as coordinates of each point that makes up the polygon, the boundary box in which the polygon will be plotted and the Coordinate Reference System (CRS) that is used to calculate where the polygons should be placed on the map.

As such, we need to bring in another file that contains this data.

So, lets read in a shapefile (.shp) that contains the polygons we need. Because shapefiles are more complex than simple tables they need to be loaded into an object of class ‘SpatialPolygonsDataFrame’. We do this by loading the file in using a specific function: ‘readOGR’. The ‘readOGR’ is a function of the ‘rgdal’ package.

The complex nature of the SPatialPolygonsDataFrame can be explored by looking at the different parts that make them up. Each SpatialPolygonDataFrame has ‘slots’ that each contain differnt parts. A basic way to explain their structure can be to think of them as a datafram of lists, with each list being different to the others.

Before doing this you will need to ensure a copy of the shapefile is in the working directory as the ‘readOGR’ function requires a reference to the file location.

Lets load it into an object (‘shapes’) and take a look at the names of the ‘slots’;

shapes <- readOGR(dsn=getwd(), layer = "Regions_December_2018_EN_BFC")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Anthony\Documents\GitHub\Leaflet_examples", layer: "Regions_December_2018_EN_BFC"
## with 9 features
## It has 9 fields
## Integer64 fields read as strings:  objectid bng_e bng_n

slotNames(shapes)

## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

Notice that there are the following ‘slots’; - data: data related to the polygons - polygons: point coordinates of each point that makes up the polygon - plotorder: a simple mnemonic that dictates the order in which the polygons will be plotted - bbox: the boundary box in which the polygon will fit - proj4string: information on the CRS and projection of the polygons

Let’s take a look at the structure of the SpatialPolygonsDataFrame as this will help you understand the next steps we need to take before the map can be plotted.

glimpse(shapes)

## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  9 obs. of  9 variables:
##   .. ..$ objectid  : Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9
##   .. ..$ rgn18cd   : Factor w/ 9 levels "E12000001","E12000002",..: 1 2 3 4 5 6 7 8 9
##   .. ..$ rgn18nm   : Factor w/ 9 levels "East Midlands",..: 4 5 9 1 8 2 3 6 7
##   .. ..$ bng_e     : Factor w/ 9 levels "285015","350015",..: 4 2 5 7 3 9 8 6 1
##   .. ..$ bng_n     : Factor w/ 9 levels "102567","172924",..: 9 8 7 6 5 4 3 2 1
##   .. ..$ long      : num [1:9] -1.73 -2.77 -1.29 -0.85 -2.2 ...
##   .. ..$ lat       : num [1:9] 55.3 54.4 53.9 52.8 52.6 ...
##   .. ..$ st_areasha: num [1:9] 8.59e+09 1.42e+10 1.54e+10 1.56e+10 1.30e+10 ...
##   .. ..$ st_lengths: num [1:9] 1008416 2075325 1385385 1332316 953902 ...
##   ..@ polygons   :List of 9
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   ..@ plotOrder  : int [1:9] 9 6 8 4 3 2 5 1 7
##   ..@ bbox       : num [1:2, 1:2] 82672 5338 655645 657534
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

glimpse(shapes@proj4string)

## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps"| __truncated__

The shapefile is projected in a different CRS to the standard CRS used by leaflet: transverse mercator (tmerc) using the GB36 datum. This information is stored in this section;

proj4string :[+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs

By default, OS boundaries are stored in the GB1936 CRS. We can either change the CRS being used by leaflet or transform the polygons data to a more faviourable CRS. Transforming is easier, lets change the CRS on the polygons to the same as the standard CRS in leaflet - WGS84 (also known as EPSG:4326).

To transform one CRS to another we can use the ‘spTransform’ function of the ‘sp’ package and we’ll check that the transformation has been completed;

shapesWGS84 <- spTransform(shapes, CRS('+init=epsg:4326'))
glimpse(shapesWGS84@proj4string)

## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

plot(shapesWGS84)

So, we have our polygons ready to plot on a map but we are still not ready to take this step just yet. Every shapefile has a data portion that contains data about the polygons. The content of the data portion varies depending on the shapefile you are using. We can take take a closer look the structure the data portion by referencing that slot with ‘@’;

str(shapesWGS84@data)

## 'data.frame':    9 obs. of  9 variables:
##  $ objectid  : Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9
##  $ rgn18cd   : Factor w/ 9 levels "E12000001","E12000002",..: 1 2 3 4 5 6 7 8 9
##  $ rgn18nm   : Factor w/ 9 levels "East Midlands",..: 4 5 9 1 8 2 3 6 7
##  $ bng_e     : Factor w/ 9 levels "285015","350015",..: 4 2 5 7 3 9 8 6 1
##  $ bng_n     : Factor w/ 9 levels "102567","172924",..: 9 8 7 6 5 4 3 2 1
##  $ long      : num  -1.73 -2.77 -1.29 -0.85 -2.2 ...
##  $ lat       : num  55.3 54.4 53.9 52.8 52.6 ...
##  $ st_areasha: num  8.59e+09 1.42e+10 1.54e+10 1.56e+10 1.30e+10 ...
##  $ st_lengths: num  1008416 2075325 1385385 1332316 953902 ...

We can join the accident data we prepared to our polygons data using a simple join. But in order to do that, we are going to need to clean up the data portion first;

# region was read in as a factor, which will affect any joins.  change region data type from factor to char
shapesWGS84@data$rgn18cd <- as.character(shapes@data$rgn18cd)

# also, lets rename the region column for ease of use
colnames(shapesWGS84@data)[2] = "rgn"

Now, we can do our join! The follwoing join will bring in all columns from our accident data, although in this case there is only the ‘count’ column to bring;

shapesWGS84@data <- shapesWGS84@data %>%
  left_join(uk_accidents_2017, by = "rgn")

# check it worked
head(shapesWGS84@data)

##   objectid       rgn                  rgn18nm  bng_e  bng_n      long
## 1        1 E12000001               North East 417313 600358 -1.728900
## 2        2 E12000002               North West 350015 506280 -2.772370
## 3        3 E12000003 Yorkshire and The Humber 446903 448736 -1.287120
## 4        4 E12000004            East Midlands 477660 322635 -0.849670
## 5        5 E12000005            West Midlands 386294 295477 -2.203580
## 6        6 E12000006          East of England 571074 263229  0.504146
##        lat  st_areasha st_lengths count
## 1 55.29703  8592429312  1008415.7  4088
## 2 54.44945 14164090941  2075324.9 13164
## 3 53.93264 15409123403  1385384.9 11336
## 4 52.79572 15643271928  1332315.8  8734
## 5 52.55697 13003737767   953901.5 10977
## 6 52.24067 19135563220  2646523.2 12484

Success! Our ‘count’ column is at the very end of the data! Now, if we plot our polygons, the ‘count’ data will be available for us to reference on the map.

A good example of using our newly joined count of accidents in each location data is to colour each polygon based on this value. To do this we will create a colour palette using ‘colorBin’ and then add it to the polygons when they are plotted on our leaflet map using ‘addPolygons’.

# create colour palette for gradiated choropleth
accPal <- colorBin("Reds", domain = shapesWGS84@data$count)

# start adding options to customise the map;
leaflet() %>%
  addTiles() %>%
  addPolygons(data= shapesWGS84 ,weight = 1, color = ~accPal(count),
              # a label will appear when we hover over each polygon
              label = shapesWGS84@data$rgn18nm,
              # a popup will only appear when we click on a polygon
              popup = ~paste0( "<h3>", rgn18nm, "</h3>", "No.of Accidents: ", count),
              # highlights are shown on a polygon when you hover over it
              highlight = highlightOptions(weight = 2, color = "blue", bringToFront = TRUE)
              )

That’s all folks! :)